/* ================================================================================================

hrn.hpp

author: Antoine Allard
affil.: Universite Laval, Quebec City, Quebec
www:    dynamica.phy.ulaval.ca


last modif. August 27, 2013

// ==============================================================================================*/

#ifndef HRN_HPP_INCLUDED
#define HRN_HPP_INCLUDED

// standard library
#include <string>
#include <vector>
#include <cmath>

#include "matrix.hpp"

using namespace std;


class degree_distribution: public matrix
{
  private:
    // average degree of nodes
    double average_degree;
    // the distribution itself
  	vector< vector< double> > degree_dist;
    // the relative weights of each coreness
    vector<double> w_c;
    // average degree of each coreness
    vector<double> k_c;
  public:
    // constructor
    degree_distribution(string);
    //
    double operator()(int c, int k){return degree_dist[c][k];}
    //
    double get_w(int c){return w_c[c];}
    //
    double get_class_average_degree(int c){return k_c[c];}
    //
    double get_average_degree(){return average_degree;}
};

degree_distribution::degree_distribution(string _mat_filename)
{
  // sets the filename of the matrix
  mat_filename = _mat_filename;
  // determines the bounds of the degree distribution
  set_dimensions();
  // copy the elements of the matrix from the file to the object
  fill_matrix(degree_dist);
  // transforms the K matrix into the degree distributions (the same object is reused to save memory)
  norm_matrix(degree_dist);
  w_c = vector<double>(dim0,0.0);
  for(int c(0); c<dim0; ++c)
  	for(int k(0); k<dim1; ++k)
  	  w_c[c] += degree_dist[c][k]; 
  norm_row_matrix(degree_dist);
  // computes the average degree of each coreness class
  k_c = vector<double>(dim1,0.0);
  for(int c(0); c<dim0; ++c)
    for(int k(0); k<dim1; ++k)
      k_c[c] += k * degree_dist[c][k];
  // computes the average degree of nodes
  average_degree = 0.0;
  for(int c(0); c<dim0; ++c)
    average_degree += w_c[c]*k_c[c];
}


class transition_probabilities: public matrix
{
  private:
    //
    vector< vector< double > > C_mat;
    //
    vector< vector< vector< vector< double > > > > transition_prob;
  public:
    // constructor
    transition_probabilities(degree_distribution&, string);
    //
    double operator()(int to_c, int to_s, int from_c, int from_s){return transition_prob[to_c][to_s][from_c][from_s];}
};

transition_probabilities::transition_probabilities(degree_distribution& DD,string _mat_filename)
{
  // sets the filename of the matrix
  mat_filename = _mat_filename;
  // gets the bounds of the matrix
  set_dimensions();
  // copy the elements of the matrix from the file to the object
  fill_matrix(C_mat);
  // adds the transpose to the matrix to make it symetric
  add_own_transpose(C_mat);
  // normalizes its elements to one
  norm_matrix(C_mat);
  // initialization of "transition_prob"
  transition_prob.resize(dim0);
  for(int c(0); c<dim0; ++c)
  {
    transition_prob[c].resize(2);
    for(int i(0); i<2; ++i)
    {
      transition_prob[c][i].resize(dim0);
      for(int cp(0); cp<dim0; ++cp)
      {
        transition_prob[c][i][cp] = vector<double>(2,0.0);
      }
    }
  }
  // builds the transition probabilities
  double sum_under, sum_over, frac_red, frac_blue;
  for(int c(1); c<dim0; ++c)
  {
    if(DD.get_w(c) > 1e-12)
    {
      // pre-calculation
      sum_under = 0.0;
      for(int cp(0); cp<c; ++cp)
        sum_under += C_mat[c][cp];
      sum_over  = 0.0;
      for(int cp(c+1); cp<dim0; ++cp)
        sum_over += C_mat[c][cp];
      frac_red  = DD.get_w(c) * c / DD.get_average_degree();
      frac_blue = DD.get_w(c) * (DD.get_class_average_degree(c) - c) / DD.get_average_degree();
    
      if(frac_red > 1e-12)
      {
        // equation (3a)
        for(int cp(c+1); cp<dim0; ++cp)
          transition_prob[cp][1][c][0] = C_mat[c][cp] / frac_red;
        // equation (3e)
        transition_prob[c][1][c][0] = (frac_blue - sum_under) / frac_red;
        // equation (3f)
        transition_prob[c][0][c][0] = (2*frac_red - C_mat[c][c] - 2*sum_over) / frac_red;
      }
      if(frac_blue > 1e-12)
      {
        // equation (3b)
        for(int cp(1); cp<c; ++cp)
          transition_prob[cp][0][c][1] = C_mat[cp][c] / frac_blue;
        // equation (3d)
        transition_prob[c][0][c][1] = (frac_blue - sum_under) / frac_blue;
      }
    }
  }

  // clears the C matrix
  C_mat.clear();
}


double g(vector< vector< double > >& x, double T, int c, degree_distribution& P, transition_probabilities& R)
{
  // gets the dimensions of the degree distribution
  vector<double> dims = P.get_dimensions();

  // evaluate the function
  vector<double> tmp_value(2,0.0);
  for(int cp(0); cp<dims[0]; ++cp)
    for(int j(0); j<2; ++j)
      for(int i(0); i<2; ++i)
        tmp_value[i] += x[cp][j] * R(cp,j,c,i);
  double return_value = 0.0;
  for(int k(c); k<dims[1]; ++k)
    return_value += P(c,k) * pow((1-T) + T * tmp_value[0],c) * pow((1-T) + T * tmp_value[1],k-c);

  return return_value;
}

double f(vector< vector< double > >& x, double T, int c, int s, degree_distribution& P, transition_probabilities& R)
{
  // gets the dimensions of the degree distribution
  vector<double> dims = P.get_dimensions();

  // evaluate the function
  vector<double> tmp_value(2,0.0);
  for(int cp(0); cp<dims[0]; ++cp)
    for(int j(0); j<2; ++j)
      for(int i(0); i<2; ++i)
        tmp_value[i] += x[cp][j] * R(cp,j,c,i);

  double return_value = 0.0;
  if(s == 0) // red stub
    if(c != 0)
      for(int k(c); k<dims[1]; ++k)
        return_value +=  P(c,k) * pow((1-T) + T * tmp_value[0],c-1) * pow((1-T) + T * tmp_value[1],k-c);
  if(s == 1) // blue stub
    if(c != 0)
    {
      for(int k(c+1); k<dims[1]; ++k) // the case k=c does not contribute and is therefore excluded to avoid 0 times infinite
        return_value += (k-c) * P(c,k) * pow((1-T) + T * tmp_value[0],c) * pow((1-T) + T * tmp_value[1],k-c-1);
      if((P.get_class_average_degree(c) - c) > 1e-12)
        return_value /= P.get_class_average_degree(c) - c;
    }

  return return_value;
}

double compute_S(double T, degree_distribution& P, transition_probabilities& R)
{
  // dimensions of the degree distribution
  vector<double> dims = P.get_dimensions();
  
  // initial conditions
  vector< vector< double > > old_a(dims[0],vector<double>(2,0.0));
  vector< vector< double > > new_a(dims[0],vector<double>(2,0.0));
  
  // finds the fixed point
  bool new_it_needed = true;
  while(new_it_needed)
  {
    // computes the next point
    for(int c(0); c<dims[0]; ++c)
      for(int i(0); i<2; ++i)
        new_a[c][i] = f(old_a,T,c,i,P,R);

    // new iteration needed?
    new_it_needed = false;
    for(int c(0); c<dims[0]; ++c)
      for(int i(0); i<2; ++i)
        if(abs(new_a[c][i]-old_a[c][i]) > 1e-8)
          new_it_needed = true;

    // make a copy of the actual position
    for(int c(0); c<dims[0]; ++c)
      for(int i(0); i<2; ++i)
        old_a[c][i] = new_a[c][i];
  }

  // computes the size of the giant component
  double S = 1.0;
  for(int c(0); c<dims[0]; ++c)
    S -= P.get_w(c) * g(old_a,T,c,P,R);

  return S;
}

#endif // HRN_HPP_INCLUDED